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Abstract 

We present an exploratory study of chiral effective theories of nuclei with methods 
adopted from lattice quantum chromodynamics (QCD). We show that the simulations 
in the Euclidean path integral approach are feasible and that we can determine the 
energy of the two nucleon state. This opens up the possibility to determine in future 
simulations nucleon phase shifts by varying the parameters and the simulated volumes. 
The physical cut-off of the theory is realised by blocking of the lattice fields. By keeping 
the block size fixed in physical units the lattice cut-off (i.e. the lattice spacing) can be 
freely changed. This offers an effective way for controlling lattice artefacts. 
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1 Introduction 



Quantum Chromodynamics (QCD) is the nowadays accepted theory of strong interactions 
in terms of the fundamental quark and gluon degrees of freedom. Also, if one is interested 
in nuclear physics, QCD is the relevant theory to calculate for instance binding energies of 
nuclei from first principles. However, QCD is strongly interacting at low energies and hence, 
non-perturbative methods are needed in order to study quantities as, for instance, nuclear 
binding energies. 

The method of choice for ab initio, non-perturbative investigations of QCD is provided 
by lattice QCD. Even if this field was facing tremendous progress over the last few years, the 
computation of nuclear binding energies for say A > 2 is still out of reach. But it is also not 
clear whether those quantities need to be computed in the fundamental theory, because the 
relevant degrees of freedom are effectively given by hadrons described in a chiral effective 
theory [HE], and not by quarks and gluons in QCD. And again, lattice methods can be 
used in order to investigate these chiral effective theories from first principles, for a review 
see ref. p]. For an overview over the various available approaches and methods to nuclear 
physics using chiral effective theories we refer to recent review articles [11|5]. 

A lattice approach to nuclear physics was presented in a row of papers by Borasay and 
collaborators P4T2]. The authors of these publications take the approach of simulating the 
chiral effective theory on a discrete space-time lattice. Starting point for these simulations 
is the non-relativistic effective theory. As a consequence. Quantum Monte-Carlo methods 
based on Fock states [IS] are used for the numerical simulations. The lattice regularisation is 
utilised as the physical cut-off in the theory, hence it cannot be removed. It naturally takes 
values of a ~ 2 fm. Lattice artifacts are reduced perturbatively by adding counter-terms 
to the simulated action. Such large values of the lattice spacings allow for simulations with 
L = 8 lattice points at most in each of the spatial directions. This is important, as in the 
simulations a sign problem is apparent, which is controllable for small L, but severe for large 
L. 

Apart from the fact that QCD and the chiral effective field theory (ChET) are rather 
different, also the approach taken in these publications differs fundamentally from the meth- 
ods applied in LQCD simulations. In LQCD the path integral in Euclidean space-time is 
directly simulated using appropriate Markov chain Monte Carlo methods. Simulations are 
performed in a finite volume with discretised space-time, with the inverse lattice spacing 
serving as a momentum cut-off. Renormalisaton is then performed by properly removing 
the cut-off in the continuum limit when the lattice spacing tends to zero (a — )■ 0). At this 
point the fundamental difference between QCD and ChET becomes evident: LQCD has a 
non-trivial continuum limit and hence the cut-off can be completely removed. In ChET the 
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cut-off cannot be removed - in fact it has a physical relevance representing the energy scale 
where the basic degrees of freedom (pions and nucleons) become inappropriate for describ- 
ing physics. In a perturbative framework the necessity of a finite cut-off is implied by the 
non-renormalisability of the theory. In a non-perturbative formulation ChET belongs to the 
general class of Yukawa-type theories with fermion-boson interactions. (In a four-fermion 
theory, for instance, the boson in the Yukawa interaction is generated dynamically.) Yukawa- 
theories are expected to have only a trivial (non- interacting) continuum limit, therefore the 
cut-off is necessary for maintaining a non-zero interaction. (For an introduction in lattice 
Yukawa models see Chapter 6 in [H].) The finite value of the cut-off also corresponds to 
the non-zero size of pions and nucleons below which distance scale obviously some other 
description is required (namely, with quarks and gluons in terms of QCD). 

In the continuum limit of LQCD all symmetries of the continuum theory broken by 
discretisation, including Lorentz symmetry, are restored and lattice artifacts in all quantities 
are removed. Performing the continuum limit requires an extrapolation to vanishing value 
of the lattice spacing a — )■ and therefore the simulations have to take place in ranges of 
a which allow such an extrapolation in a controlled way. It is not a priori clear how small 
values of the lattice spacing are required for this procedure, and eventually one will only find 
it out empirically. 

The so called Symanzik effective theory p^HTT] allows for a better understanding which 
lattice artifacts one has to expect. For sufficiently small values of a LQCD can be described 
by an effective local action [TS] 



where 5*0 is the continuum QCD action and the additional terms Sk for A; > are to be 
interpreted as operator insertions in the continuum theory 



with Ck being combinations of local composite fields with mass dimension A + k. The list 
of possible fields is constrained by the symmetries. Similarly one can write for a local gauge 
invariant composite field 0(x), where we neglect mixing under renormalisation for simplicity. 



The fields 0^ are again linear combinations of local fields with appropriate symmetries and 
dimension. Without discussing all the details - which can be found in ref. [IS] - one finds 



S'eff — 5*0 + aSi -|- a 5*2 + . . . , 



(1) 





(2) 
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that a connected n-point function on the lattice to leading order in the lattice spacing reads 

-a j (0o(xo)---0o(x„)/:i(y))^°'^* 

(3) 



+ a ^(0o(a;i) ■ ■ ■ 0i(xfc) ■ ■ ■ 0o(a;„))'°"* 
fc=i 

+ ... . 



From this expression it should be clear that first of all the value of a should be small enough 
to neglect all higher order terms in the lattice spacing. In particular, on the right hand 
side of eq. ([3]) there should be only combinations with a{. . .) <^ 1. This is also the reason 
why LQCD investigations with heavy quarks appear to be difficult since then arriq can easily 
become of order 1 and can lead to significant lattice artifacts. Of course, all terms on the 
right hand side of eq. Q come with unknown coefficients, so also with aniq ^ 1 large lattice 
artifacts are not automatically to be expected. 

The so called Symanzik expansion can also be used to systematically remove lattice 
artifacts from the theory. This requires to add counter terms to the lattice action with 
coefficients that need to be determined non-perturbatively. Improving the action is, however, 
not always enough. There are observables that require operator specific improvement with 
coefficients that again need to be determined non-perturbatively. 

Comparing the lattice QCD approach and the lattice approach to chiral effective theories 
from refs. [6HT2] rises an immediate question, namely how large are lattice artifacts in the 
lattice simulations of the chiral effective theory and whether they are will controlled. This 
question is hard to answer since the continuum limit a — )■ cannot be performed due to the 
physical interpretation of the cut-off. Moreover, a variation of the cut-off is for technical 
reasons difficult, but a variation in a reasonable range is desirable. 

This is why we take a different approach to lattice simulations of these chiral effective 
theories, which we shall present in this paper. It consists of using the fully relativistic path 
integral formulation of the chiral effective theory in Euclidean space-time, which can be 
simulated by means of Markov chain Monte Carlo methods. The physical cut-off of the 
theory is implemented by block-field methods. The sizes of the blocks represent the hadron 
sizes. These have to be kept fixed if the lattice spacing is changed, for instance, keeping 
R ■ M fixed where R denotes the block size and M is a physical mass, say the nucleon 
mass. Separating the lattice cut-off from the physical cut-off allows to change the lattice 
cut-off (i.e. the lattice spacing) by keeping the physical content of the theory unchanged. On 
the basis of the Symanzik effective theory it has to be expected that by making the lattice 
spacing sufficiently small the lattice artifacts, for instance, Lorentz-symmetry breaking can 



be suppressed. One can also imagine to reach in the hmit a — )■ a non-trivial continuum limit 
describing an interacting non-local quantum field theory. Strictly speaking, the existence of 
such a continuum limit is not known at present - similarly to the existence of the continuum 
limit of many other lattice quantum field theories. 

A substantial simplification in fermionic lattice quantum field theories is the so called 
quenched approximation which corresponds to omitting closed fermion loops. Since both 
neutron and proton are heavy, we do not expect that nucleon-anti-nucleon loops play a 
significant role and hence perform the simulations in the quenched approximation. Later on 
one can, of course, perform simulations in the full theory, removing this approximation. 

Another possible simplification is to work within a non-relativistic approximation - in- 
stead of in a relativistic quantum field theory as we are doing here. One big advantage of 
the non-relativistic formulation might be that the rest masses of the nuclei do not appear in 
the theory and binding energies are computed directly, whereas in our relativistic approach 
binding energies are sub-percent effects as compared to the energy levels we determine. This 
is one point we shall address in a forthcoming publication. 

The authors of refs. P^VET] follow a similar approach to the one presented here, but 
applied to a Yukawa model with one scalar field. Previous investigations of Yukawa models 
in the Euclidean path integral formulation have been performed in the context of upper and 
lower limits of masses in the electroweak theory, see for instance refs. [22|l23]. 

In this paper we shall discuss our method and show evidence for its applicability. Physical 
results are planned to be presented in forthcoming publications. In the following section we 
discuss our lattice action followed by details of the numerical methods we apply. We shall 
close with presenting some simulation results and give a conclusion and outlook. 

2 Lattice actions 

Let us start by introducing our notations. The nucleon fields are described by a pair of 
Grassmann variables and are denoted by 

'4'ax , i'ax (a = 1, 2) (4) 

where a is the isospin doublet index and x denotes the lattice sites. The Dirac-index of 
the nucleon field is not displayed here and the isospin index of the nucleon field will also be 
omitted in most formulae. The real boson fields are 

(a = 1,2, 3) (5) 

with the triplet isospin index a. Here cp^'^^ stands for the pion and 0^°-' and (p^^^ are Hubbard- 
Stratanovich auxiliary boson fields with isospin zero and one, respectively. 
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These latter are used to describe four-nucleon interactions which in the Euchdean lattice 
action have the following form: 



^NA = Yl + + Co 0f iiJx^Px) + C, iij^TaiJx) } , (6) 

where Tq, (a = 1, 2, 3) are Pauli-matrices for isospin and a summation over repeated indices a 
is understood. The four-nucleon interactions are obtained after integrating over the auxiliary 
fields according to 

# exp{-02 - = exp{^ . (7) 

Our choice of the pion-nucleon interaction corresponds to the lattice discretisation of the 
continuum interaction d^iiai.'^^lbln'Ta'^) and is given as 

= E { T ('^-+'^ ~ '^--a) (^^^57.r„^.) I , (8) 

X ^ ■' 

where fi denotes, as usual, the unit vector on the lattice in direction ^{fi = 1,2,3,4). 
(Similarly to isospin index a, over repeated direction index /i a summation is implied.) The 
above expression corresponds to a particularly simple discretisation of the derivative of the 
pion field which can, of course, be chosen differently. 

The parameters of the lattice action are always dimensionless. The connection to the 
(eventually) dimensionful parameters in continuum formulations is established with the mul- 
tiplication by an appropriate power of a mass parameter. 

Let us illustrate this on the example of the four-nucleon interactions in ([n])-(IZl)- The 
(bare) nucleon field ipx is related to the (bare) continuum nucleon field by ipx = a^^'^ipT^^ (and 
similarly for ipx)- The relation for the auxiliary scalar fields is given by (px = 0^0^°"^*. This 
implies that the four-nucleon couplings Co,i are related to their continuum counterparts by 
Co,i = a~^CQ°^^ = {mj<fCQ"^^) / {amjsf) . Here, mNC™]^* is again dimensionless and independent 
of the lattice spacing a. In the following we always work in lattice units. The connection 
to the (eventually) dimensionful parameters in continuum formulations is established by 
multiplication with an appropriate power of a mass parameter and/or a power of the lattice 
spacing a. Note that in Ref. [7] the coefficients appear quadratically in the action. 

In addition to the interaction terms we also need a kinetic term for the nucleon which 
we take to be the Wilson fermion lattice action |24j : 

-^N = E ^ (^^^^) - '^N E (^^+a[1 + 7m]^^) } ■ (9) 

X i /i=±l J 
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Here is the hopping parameter defining the nucleon mass and the convention — —7^ 
is followed. For the pion field, besides the kinetic term, also a self-interaction term is allowed, 
hence we define the pion part of the action by 

where is the hopping parameter of the pion and A gives the strength of the self-interaction. 
The total lattice action is the sum of all the above terms, that is 

S — St<} + + S'na + 'S'ntt ■ (11) 



2.1 Block fields 

The interaction among hadrons, like nucleons and pions, is non-local as a consequence of 
the extended hadron structure implied by Quantum Chromodynamics. This non-locality 
can be approximately taken into account by introducing block fields in the lattice action. 
Block fields are sums of the the above local fields weighted by appropriate numerical factors. 
In order to decrease rotation symmetry breaking, our definition of the block fields tries to 
be as close as possible to an exponential decrease of the weight factors proportional to the 
Euclidean distance squared. Of course, periodic (or anti-periodic) boundary conditions have 
to be taken into account and therefore we define the squared distance between two points 
x,y on the lattice as 

{x,yf ^Y.\^„y,\' (12) 

where 

= min(|x^-y^|,|x^-y^-FL^|,|x^-y^-L^|) . (13) 

Here denotes the lattice extension in the direction /i. With this definition a generic block 
field is defined by 

Yl <PyOxp{-S {x,yf} . (14) 

y,{x,y)^<R'' 

The blocking parameters can depend on the type of fields, that is one can have different 
parameters for the nucleon (i?N, -S'n), for the pion (i?^, S'tt) and for the auxiliary fields {Ro,So) 
and (i?i, 5"!), respectively. The lattice action in terms of the block fields has exactly the same 
form as the above action in terms of the local fields (including the summation over the lattice 
sites). The only changes are: 

V^,^*,, {A = n,Q,l). (15) 
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Besides reflecting the non-locality of the interactions, the blocking also has an important 
role in defining the cut-off of the theory. The high momentum modes are decoupled from the 
interactions by the blocking. Only the modes below some momentum cut-off are interacting, 
the cut-off value being determined by the blocking parameters {R, S). The cut-off introduced 
by the blocking can be considered as a physical effect. The lattice cut-off can be changed 
independently of it. This allows to move the lattice spacing to small values for reducing 
lattice artefacts. 

3 Numerical simulations 

3.1 Lattice parameters 

In order to gain experience with the lattice formulation defined in the previous section we 
started numerical simulations on small lattices and also introduced some simplifications in 
the choice of lattice action parameters. In most cases we did not block the nucleon field, 
only the bosonic fields and fixed the blocking parameters as follows: 

i?N = , So = Si = S^ = 2.0, Rq = R^ = R^ = 1.5 . (16) 

This choice of the radius sets the number of points in a block to be 33. 

The simplest choice of couphng parameters is to keep only Cq and Ci defining the four- 
nucleon interactions. We also experimented with the additional introduction of the pion field 
and the coupling but here we report only on results in the pion-less theory. A detailed 
investigation of the effect of the pion coupling is left to future work. Following Ref. [7] Cq is 
taken to be real which corresponds to an attractive interaction whereas Ci = i |Ci| is purely 
imaginary describing a repulsive interaction. In most cases the magnitudes of Cq and Ci 
were equal, that is 

a = 0, Co = |Co|, Ci = z|Ci|, |Co| = |Ci|. (17) 

(Note that according to eq. ([7]) the signs of Cq and Ci are irrelevant.) 

For this exploratory work we had access to a PC cluster with 12 graphics processing units 
(CPUs) attached to it. Due to the memory limitations of the CPUs the maximally feasible 
spatial lattice size was 32^, {L = 32). The temporal lattice extension is taken four times 
longer in order to allow for a precise determination of the energy values. In summary, we 
did simulations on three types of lattices: 

16^-64, 24^-96, 32^-128. (18) 
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This restriction to the currently feasible lattice volumes has implications on the choice of 
parameters we can simulate: in order to study the two nucleon system a large physical 
volume is required, because otherwise the smallest available lattice momentum is too large. 
On the other hand we would like to simulate with as small as possible values of the lattice 
spacing a. To set the scale we fixed the equal mass of the proton and neutron to be Mn = 
-^proton = ^neutron = 939MeV/c^. This means that, for instance, a nucleon mass in lattice 
units oMn = 1 implies a lattice spacing a ~ 0.21 fm and in our cases spatial extensions of 
16 a ~ 3.4 fm, 24 a ~ 5.0 fm and 32 a ~ 6.7 fm, respectively. In order to reach larger volumes 
one can, of course, increase the nucleon mass in lattice units but in this way lattice artefacts 
will also increase. 

With these choices of volumes and with aM^ = 1 we get for the smallest non-zero 
momentum (27r)/(16La) ~ 369MeV/c, (27r)/(24La) ~ 246MeV/c and (27r)/(32La) ~ 
184MeV/c, respectively. For extracting physics using the approach presented here both 
lattice sizes and minimal non-zero momenta would require a factor, say, four increase of the 
lattice extensions which will be comfortably possible to reach with present day computer 
resources. 

3.2 Numerical methods 

Our present aim is to determine the energies (masses) of different nuclear systems. This 
can be achieved by investigating the large (Euclidean-) time behaviour of different sets of 
correlators. Since at present we restrict ourselves to quenched simulations, where the fermion 
determinant of the nucleon field is neglected, the creation of the boson field configurations is 
simple. In case of the auxiliary fields and $1:^^ one has to create Gaussian distributions. 
The pion field can also be simply produced by some update algorithm as, for instance. 
Metropolis algorithm - the only mild complication being to take into account the non-locality 
introduced by the blocking. 

The nucleon mass can be determined from the behaviour of the nucleon time-slice cor- 
relator. The time-slice operators are defined as 

X4=t X4=t 

= ^ ^ ^xi ,a;2 jXs ,a;4 ) -^t = ^ ^ ^xi ,a;2 ,0:3 ,a;4 (l^) 

Xj jX2 jX'J, X\ jX2 ;^'3 

and the nucleon correlator is, with a Dirac-projection to the state propagating in positive 
direction, 

TrDirac[(l+74)(iVtliVi,)] . (20) 

The expectation value of the fermion bilinear gives a fermion propagator which is the inverse 
of the fermion matrix in the fermionic part of the action. The overwhelming part of computer 

9 



resources in our quenched simulations is spent in the calculation of the fermion propagators 
by an iterative inverter of this sparse matrix. 

The computation of nucleon propagators has been done most of the time by applying a 
mixed precision Conjugate Gradient inverter, see the appendix. The crucial problem for the 
inverter is to deal with the very small values of the nucleon propagators at large distances. 
The solution of this problem is to use distance preconditioning following Ref. [25]. Since the 
nucleon propagator behaves nearly exponentially for distances which we use for extracting 
the masses (in most cases up to a time distance half the time extension L4 of the lattice), 
we choose the preconditioning function to be 

( exp{-Pt} if t < L4/2, 

a{t) = \ (21) 
[ exp{-P(L4-t)} if t>U/2. 

The parameter P can be chosen typically by an amount 0.1 — 0.5 smaller than the nucleon 
mass in lattice units am^. 

In order to obtain the masses of multi-nucleon (in the present paper two-nucleon) states 
with sufficient precision, one has to find the proper composite operators defining the correl- 
ators. Here we restrict ourselves to proton-neutron states. For local operators we take in 
the spin-0 and spin-1 channels, respectively, 

^ixC75^2x, *ixC7fc^2x, (A; = 1,2, 3), (22) 

where C denotes the charge conjugation Dirac matrix. 

Especially for scattering states it is important to also take extended (smeared) operators 
where the proton and neutron are at different points. In case of Gaussian smearing one can 
use the smearing function 

exp {-cri|xi,yip - o-2k2,l/2p - CTsksjI/sP} (23) 

with the notation introduced in ( IT3|) . For a spherical state in the spin-0 channel one can 
take a = 0"! = (T2 = a^. For spin-1, ellipsoidal states with e.g. Ci ^ a2 = cr^ are useful. In 
this latter case we also tried linear smearing corresponding to cr2,cT3 = 00. In order to save 
computer time one can cut the summation over sites off at distances where the smearing 
function in (1231) is smaller than, say, 10~^. In case of spherical smearing this corresponds to 
a cut-off radius of 

p={MM}"\ (24) 

The simplest way to determine the masses is to fit some of the correlators by an expo- 
nential function in time intervals for distant time-slices. In case of small enough statistical 
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errors one can also obtain good fits with a sum of two (or more) exponentials. The best 
results can be achieved, however, by taking a set of some operators in a given channel and 
calculate the correlator matrix among them. For determining the energies of two-nucleon 
(actually proton-neutron) states we typically start from a 4 x 4 correlator matrix. The four 
states are chosen from local, spherically smeared and elliptically smeared states with different 
Dirac-matrices. 

The correlator matrix can be approximated by the sum of contributions of eigenstates of 
the Hamiltonian (i. e. of the transfer matrix). In general, a real symmetric D x D correlator 
matrix C{t2,ti) between time-slices ti and t2 > ti is defined by the matrix elements of 
D operators Oa, Ob, ■ ■ ■ , Od- If the energy eigenstates are |n), n = 1,2, ... ,M then in a 
shorthand notation 

^ C{t2,ti)aa C{t2,ti)ab ■■■ C'(t2,^l)ad^ 
, x_ C{t2,ti)ab C{t2,ti)bb ... C{t2,ti)bd 

\ C(t2,ti)ad C(t2,ti)bd ••• C(t2,ti)dd J 

where the matrix elements can be written as, for instance, 

Cit2,ti)ab = ia\l)t,m)t, + ia\2)t,m)t, + ... + ia\M)t,{b\M)t, (26) 

with 

{c\k)t = {0\OMk) = {k\0,{tm , (27) 

for c = a,b, . . . ,d and k = 1,2, . . . , M. 

Assuming that we consider bosonic (fermionic) operators, we have periodic (anti-periodic) 
time dependence with the time extension of the lattice L4. This implies 

{a\k)t,{b\k)H = ia\k) ib\k) {exp[-tEk] ± exp[-(L4 - t)Ek]} . (28) 

where the positive and negative sign stands for periodicity and anti-periodicity, respectively. 
Here t = t2 — ti, is the energy (e.g. mass) corresponding to the state \k) and 
[a\k) = (a|/c)o, {b\k) = (&|fc)o- Fitting the correlator matrix by the expression given by fl2S]) 
- (128|) one can obtain the energies we are looking for |26] . The statistical errors of the results 
can also be obtained by methods similar to those described in Section 5 of this reference. 

Since in the present case the relevant (multi-) nucleon correlators can be determined to 
a very good precision, one can perform least-square fits by minimising the correlated chi- 
squared. In order to obtain a good starting point for the minimisation, one can first minimise 
the uncorrelated chi-squared defined by 



An 



E(^^) (29) 



i=l 



5X, 
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where the index i runs over the independent matrix elements to be fitted, Xi and 5Xi are the 
mean value and error of the matrix element respectively, and fi{p) is the fitting function 
of Np parameters {pi,p2, . . . ,Pnp) defined by The best fit obtained in this way 

can be taken as a starting point to minimise the correlated chi-squared 

Nc 

Xl = E (/^(P) - ^0 {fj{p) - X,) , (30) 

where Mij = NC~^ ^ with the number of input data and the correlator matrix 

1 ^ — — 

n=l 

In general, the correlator matrix in ( l3T]) can be determined with sufficient precision for 
obtaining its inverse and its eigenvectors. In some cases, in particular if the dimension of the 
correlator matrix Nq is large, smoothing of the smallest eigenvalues j^ZlEH] can be helpful 
but does not substantially influence the results. The advantage of properly obtaining the 
minimum of xl is that one can select "good fits" by the value of xl per number of degrees 
of freedom {Nc — Np). The mean value and error of a quantity is defined by considering the 
distribution of its values in good fits. The quoted value is then the position of the median 
of the distribution of these selected values. The error defines a (symmetric) interval around 
the median such that 68% of the distribution is contained in it. 



3.3 Numerical results 

The physical quantities we are interested in are for instance the nucleon-nucleon scattering 
length and binding energies of multi nucleon states. In a lattice simulation, the determ- 
ination of these quantities requires a study of the (finite) volume dependence of one and 
two (and multiple) particle energies. In this methodical paper we hence try to understand 
how precisely the corresponding quantities, i.e. the nucleon and two-nucleon masses, can be 
determined. 

In order to do so, we performed several simulations and determined the masses as de- 
scribed in the previous sub-section. A typical example is a run on a 32^ ■ 128 lattice at 
kn = 0.08, Co = 0.2, Ci = 0.2i. It turned out that the masses can be very precisely 
obtained even from a modest statistical sample of 120 configurations: see figures [T]|3J 

In figure [1] we plot the actual nucleon and two nucleon correlators as functions of the time 
t in lattice units on a logarithmic scale. The decay is nearly exponential in the whole time 
range. Fitting the correlators in different time-slice distance intervals [^1,^2] by minimising 
the correlated chi-squared one finds that for the nucleon mass am-N we observe a plateau 
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Figure 1: The nucleon correlator (a) and two-nucleon correlator (b) on a logarithmic scale 
as a function of time-slice distance t in lattice units on a 32^ ■ 128 lattice at = 0.08, Cq = 
-iCi = 0.2. 

from ti = 18 on, almost independent of the value of t2 (see fig. [2]). For the two nucleon mass 
amNN the plateau sets in somewhat later (see fig. [3]). Both quantities have in common that 
the statistical errors of the single points are in the sub-percent region. The mass value and 
its error as determined from the distribution of the fit results - as described earlier - are 
indicated in both plots by the horizontal lines. 

The required precision of two-nucleon energies can be exemplified by Liischer's formula 
for extracting scattering lengths from finite size effects of the two-particle energies To 
leading order one can express the scattering length as 

aom^ = 2 . (32) 

V mN y 47r 

The masses from Figures [2ll3l namely arriN = 1.06567(25) and amNN = 2.1511(20) give for 
the right hand side a value —59(7). The physical value of the right hand side is aom-^ = 113.1. 
This corresponds to the value of the scattering length in the ^5*0 channel Cq = -1-23.76 fm and 
rriN = 939 MeV. (Here we use the sign convention for the scattering length of Ref. [2^R) 
Our value has the right order of magnitude but an opposite sign which is due to the strong 
repulsion implied by the imaginary value Ci = 0.2i. Obviously, tuning to the physical value 



^We thank the referee of our paper for drawing our attention to the different sign conventions of the 
scattering length in the hterature which we overlooked. 
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Figure 2: Results of nucleon mass am^ fits on different fitting intervals for a 32'^ ■ 128 
lattice at = 0.08, Cq = —iCi = 0.2. The fit interval [^1,^2] is specified on the x-axis by 
ti + 0.1(^2 — ti). The horizontal lines indicate the final result and error obtained from the 
distribution of the fit results. 



would require either a much smaller imaginary value or even a real value of Ci. Besides 
of this, one cannot assume that the asymptotic formula works well already on a volume of 
extension ^ 7fm. In fact, on lattices 16'^ ■ 64 and 24^ ■ 96 in the same point we obtained for 
the right hand side of f l32|) values of about 5 and 40, respectively. This shows that, in any 
case, for the determination of scattering lengths simulations on larger volumes are required. 



4 Outlook 

In this paper we have defined a theory of nucleon and pion fields in the Euclidean path 
integral formulation on the lattice. The inherent physical cut-off of this theory has been 
implemented by formulating the lattice action in terms of blocked fields. This physical cut- 
off is given by the block size. In this way the lattice cut-off is separated from the physical 
cut-off and can be changed in order to determine the size of lattice artefacts. 

The positive outcome of the first studies we have performed is that the correlators of 
single- and two-nucleon systems can be determined very precisely in order to obtain the 
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Figure 3: The same as FiglEfor two-nucleon masses amNN- 

energies to a very good precision. In fact a few thousands of configurations are sufficient for 
a precision better than one per mill. Using Liischer's formula we could determine values for 
the nucleon-nucleon scattering length. The corresponding results are positive in the sense 
that with larger volumes it seems to be realistic to tune its value to the physical one. 

An important step towards obtaining physical results is hence to increase the physical 
lattice sizes. As discussed in Section l3TT| a spatial lattice extension of about L = 100 would 
correspond for a nucleon mass in lattice units oMn ~ 1 to a lattice size of about La ^ 20 fm 
and a minimal lattice momentum 2it /{lOOLa) ~ 60MeV/c. 

It will be also important to introduce the pion field besides the auxiliary fields describing 
four-nucleon couphngs. 

In order to complete this sort of lattice studies an important final step is to investigate 
the dependence of the results on the lattice spacing. For this the physical parameters as, 
for instance, the block size parameters Rm^, Sm^^ and the lattice volume Lam^ have to be 
kept fixed. The couplings in the lattice action (Co,i etc.) have to be tuned for each value of 
the lattice spacing in such a way that some well chosen physical quantities (as for instance 
some nucleon phase shifts) take their physical values. Of course, if the lattice spacing gets 
smaller the required number of lattice points have to be increased correspondingly and this 
implies an increase in the required computational power. 
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A possible source of difficulties in the quenched approximation, as observed in Yukawa 
models by the authors of [T9H2T] . is the appearance of exceptional configurations with ex- 
tremely small eigenvalues of the fermion lattice action. These configurations make the de- 
termination of correlators and therefore masses practically impossible. In our case we found 
exceptional configurations for bare couplings in the range |Co|, |Ci|, \Ct,\ > 0.3. Since this 
problem does not appear in numerical simulations in Yukawa models with dynamical fermi- 
ons [2^125] , we expect that it does also disappear in our nuclear Yukawa models if dynamical 
nucleons are included in the simulation update. For real values of the couplings Cq, Ci, C-„ 
the fermion determinant is real (non-negative) therefore the known Hybrid Monte Carlo 
methods [30] can be applied in a straightforward manner. For non-real (e.g. imaginary) 
couplings the determinant becomes complex and the numerical simulation turns non-trivial, 
if not impossible. 
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A Implementation Details 

As mentioned in the introduction, we have used graphics processing units (CPUs) in order 
to perform the numerical inversions of the Dirac operator. We have 12 NVIDIA Tesla C1060 
CPUs available with four Gb of memory each. We have used NVIDIAs CUDA environment 
to implement the Dirac operator deduced from eqs. ([6]), ([8]) and Q for CPUs, which is very 
similar to available implementations for lattice QCD, see for instance ref. [31]. 

We employ a mixed precision solver using both, the CPU and the GPU. On the GPU we 
have implemented a conjugate gradient (CG) solver inverting the squared hermitian Dirac 
operator (since Ci is purely imaginary) 

The desired result is then obtained by multiplying with Q'^ . The CG solver on the GPU is 
implemented solely in single precision (32 Bit). The CG solver is called from an outer solver, 
which is run on the CPU in double precision (64 Bit). We use iterative refinement as the 
outer solver in order to solve 

Drj = (f) 
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Algorithm 1 Iterative Refinement 
Require: 0, rj, eo > 0, ej > 

1: = 

2: Tk = (j) — Drj 

3: while 110 — DxkW > eo do 

4: solve Dpk+1 = Tk for pk+i on GPU to relative precision 

5: Xfc+i = Xk+ Pk+1 

6: r^+i = - Dxk+i 
7: A; = A; + 1 
8: end while 
9: return (p 



for rj, given some source spinor field (p. The algorithm is summarised in algorithm [T] De- 
pending on the parameters kn,Cq and Ci one has to tune the precision for which to solve 
on the GPU. The usage of distance preconditioning also had significant influence on this 
tuning: the closer the preconditioning mass to the measured mass the less stable the pre- 
scribed mixed precision solver turned out to behave. Most probably due to accumulation of 
round off errors we had to reduce the number of inner iterations further and further with 
preconditioning mass approaching the measured mass. 
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